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ABSTRACT 

This  study  consists  of  the  development  of  a  computer  program  to 
numerically  solve  the  space  and  energy  dependent  multigroup  neutron 
diffusion  equations  in  a  bare  homogeneous  fast  reactor  core  or  reactor 
material  assembly. 

The  resulting  program  is  unique  in  that  it  was  designed  for  future 
use  by  Naval  Postgraduate  School  students  undertaking  experimental 
studies  in  neutron  diffusion  with  limited  time  to  determine  numerical 
solutions  for  verification  of  their  results . 

The  equations  are  solved  iteratively  in  cylindrical  geometry  using 

a  point  successive  overrelaxation  technique.     Convergence  between 

-6 
successive  iterations  was  less  than  10       after  fifty  iterations. 

The  program  was  tested  using  ANL  three  group  data.     Flux  shapes 
and  energy  spectra  were  determined  for  a  typical  fast  reactor  core 
and  for  a  solid  iron  cylinder  with  a  source  at  its  center.     The  program 
was  also  used  to  determine  criticality. 

Computation  times  were  from  one  to  ten  minutes  with  less  than 
15 OK  words  of  core  storage  using  the  IBM  360/67. 
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I.     INTRODUCTION 

A.      BACKGROUND 

It  has  been  found  that  the  energy-dependent  diffusion  equation 
provides  a  fairly  accurate  description  of  the  neutron  behavior  in  fast 
reactor  systems  [1J. 

One  way  of  conveniently  using  this  equation  is  the  now  familiar 
multigroup  method  which  was  developed  for  reactor  calculations  by 
F.  L.  Friedman,  A.   M.  Weinberg,  and  J.  A.  Wheeler  [  2]  about  1947. 
Since  that  time,  many  solutions  to  the  age-diffusion  equation  have 
resulted  from  the  use  of  the  multigroup  method.     Some  analytical 
and  semianalytical  solutions  have  been  carried  out  for  a  few  group 
equations  [3J   but  most  of  the  work  done  has  been  numerical  due  mainly 
to  the  need  for  large  numbers  of  energy  groups,  spatial  points,  and 
spatial  regions . 

The  numerical  solutions  to  unique  problems  and  classes  of 
problems  have  evolved  in  the  form  of  computer  programs  or  codes  . 
The  degree  of  complexity  and  sophistication  of  these  programs  has 
roughly  paralleled  that  of  the  digital  computer.     At  this  point  in  time, 
the  literature  abounds  with  numerous  and  varied  codes  and  references 
thereto.     Sangren  [4]    states  that  by  1959  there  existed  in. the  United 
States  more  than  300  nuclear  reactor  codes  for  digital  computers. 
Various  groups  have  classified  and  categorized  these  codes  [4]  .     In 
1961  the  American  Nuclear  Society  and  the  Argonne  National  Laboratory 
established  a  code  center  for  the  purpose  of  collecting  and  disseminat- 
ing the  multitude  of  codes  that  existed  and  for  the  future  codes  which 
were  sure  to  be  forthcoming    [5].    Around  19  65  the  European  counter- 
part to  this  center  was  established  £5]  . 

Most  of  the  more  recent  codes  (those  published  after  19  65) 
examined  in  connection  with  this  study  are  somewhat  specialized 


along  the  following  lines:    a)  number  of  energy  (lethargy)  groups; 
b)  number  of  spatial  dimensions;  c)  number  of  regions;  d)  cross-section 
data  input;  e)  geometry;  f)  number  of  mesh  points;  g)  method  of  equation 
solution;  h)  finite  difference  approximation  and  i)  many  other  features 
related  to  poisons,  burn-up,  perturbations,  computer  system,  and  others 

B.      OBJECTIVES 

This  study  was  undertaken  to  provide  a  relatively  simple-to-apply 
multigroup  diffusion  program  for  application  to  fast  reactor  cores  and 
materials.    In  the  future  students  at  the  Naval  Postgraduate  School 
will  be  able  to  use  this  program  in  conjunction  with  experimental 
neutron  diffusion  studies.     In  view  of  the  large  number  of  complex 
codes  that  exist  and  the  difficulties  related  to  modifying  these  codes 
for  use  on  a  specific  problem  and  a  specific  computer  system,  it  is 
anticipated  that  valuable  research  time  may  be  saved  by  use  of  this 
program . 

Once  the  program  was  developed,  initial  testing  would  be  done 
using  a  subcritical  assembly.     The  first  phase  of  this  test  would 
involve  the  solution  of  the  three  group  diffusion  equation  in  a  typical 
fast  core  composition.     From  this,  the  energy  and  spatial  dependence 
of  the  flux  could  be  determined  and  comparisons  made  with  the  solution 
to  the  diffusion  problem  wherein  space  and  energy  dependence  are 
treated  as  separable  variables. 

The  second  phase  of  the  initial  test  would  involve  the  problem 
of  the  typical  core  assembly  with  a  Pu-Be  neutron  source  located  at 
its  center.     This  would  then  be  modified  to  a  solid  iron  cylinder  with 
a  centrally  located  Pu-Be  neutron  source.     The  flux  shapes  and 
energy  spectra  resulting  from  this  study  could  then  be  used  for 
comparison  with  experimental  data  resulting  from  local  research. 
Experiments  in  this  area  are  currently  underway  with  applications 


to  fast  reactor  cross    section  data  evaluation,  materials  and  shielding 
reference  [  6]  . 

Using  the  program  to  determine  a  critical  composition  and  size  for 
the  typical  fast  reactor  was  another  desired  goal.    In  this  area,  a 
study  was  to  be  made  on  the  effect  mesh  spacing  has  on  criticality 
if  any. 

Other  objectives  were  to  compare  the  energy  spectra  of  the  three 
and  the  sixteen  group  solutions  and  to  ascertain  any  difference  in 
the  flux  shapes  and  spectra  by  using  the  ANL  sixteen  group  data 
and  the  LASL  sixteen  group  data. 
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II.     THEORY 

A.      THE  BOLTZMANN  TRANSPORT  EQUATION 

The  conservation  or  balance  of  neutrons  in  a  nuclear  reactor  can 
be  completely  described  only  by  simultaneously  specifying  the  distri- 
bution in  space,  time,  energy  and  direction  of  motion  of  the  neutrons 
The  general  equation  governing  this  balance  is  called  the  Boltzmann 
transport  equation  and  may  be  written  as 


PRODUCTION  -  LEAKAGE  -  ABSORPTION  = 


3>h 

at 


where   §^-    is  the  time  rate  of  change  of  neutron  density 


at 


The  neutron  density  is  in  general  described  as 


II- 1 


number  of  neutrons  at  time  t  in  volume 
dr  about  r,  whose  speeds  lie  in  dv  about 
v  and  whose  directions  of  motion  (velocity 
vectors)  lie  in  the  differential  solid  angle 
angle  d-O-  about  direction  ^n.  . 


and  neutron  flux  is  defined  as 


$Cr}  nTj  _n_,-t  )  =  ArnCr,  at,  _fL/t  ) 


II-3 


Meghreblian  and  Holmes     1     describe  Equation  II- 1  verbally  as 


number  of  neutrons  at  speed  in  dv  about 
v  moving  in  direction  -A-  which  appear 
per  unit  time  from  source  in  dr 


+11 


number  of  neutrons  of  direction  -Tl   gained 
per  unit  time  in  dr  from  scattering 
collisions  which  scatter  the  neutron  from 
all  directions  -H-    to   -H_    and  all  speeds 
to  speed  in  dv  about  v 
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-Ill 


net  number  of  neutrons  moving  in  direction 
-O-    at  speed  in  dv  about  v  that,  are  lost 
per  unit  time  by  leakage  through  the 
boundaries  of  dr 


-IV 


number  of  neutrons  at  speed  in  dv  about 
moving  in  direction  TL    which  are  removed 
per  unit  time  from  dr  by  absorption 
collisions 


-  V 


number  of  neutrons  at  speed  in  dv  about 
moving  in  direction  sx.     which  are  removed 
per  unit  time  from  dr  by  being  scattered 
into  a  new  direction 


=VI 


change  per  unit  time  in  net  number  of 
neutrons  of  direction  -?L     in  dr 


II-4 


The  first  term  in  Equation  II-4  accounts  for  the  rate  at  which 
neutrons  with  speed  in  dv  about  v  and  with  direction  in  d_f>.  appear 
from  sources  in  dr  .     This  term  is  then  written  as 

5  (r,  or}  JL,"t )  dFdftrd-?L  .  II-5 

The  second  term  of  Equation  II— 4  gives  the  total  number  of  neutrons 
which  appear  in  the  differential  element  dfdvdJL         per  unit  time 
with  speed  in  dv  about  v  and  with  direction  d-S-      about   si-       by 
being  scattered  from  all  possible  speeds  v  and  direction  of  motion  s\! 
The  number  of  neutrons  in  drdvdJrL    which  experience  this  scattering 
is 

Z.s(ru-/)  (j>ixt  rvlsL\  t  )df  dnr'dii.'  .  II-6 
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The  probability  that  the  scattering  collision  results  in  the  neutron 
having  speed  in  dv  about  v  and  direction  d  H.     about  -O-     is 
defined  by 

£  C  or  _a  ;  rv- ',  _fL'  )    = 

probability  that  a  neutron  with  an 
initial  speed  v  and  direction  of 
motion  JrL'    ,  when  scattered,  emerges 
from  the  collision  with  a  speed  in 
dv  about  v  and  direction  of  motion 
in  d  H.     about  iv     . 


II-7 


If  Equations  II— 6  and  II— 7  are  multiplied  together  and  integrated  over 
all  possible  v'  and  _n_'     ,  the  result  is  the  total  gain  of  neutrons  in 
the  differential  element  drdvd.fi-     due  to  the  scattering  in  of  neutrons 
from  higher  speeds  and  different  directions  of  motion  written 


J^,~ 


rv-SL 


II-8 


The  third  term  of  Equation  II— 4  represents  the  net  loss  of  neutrons  from 
drdJx      due  to  leakage  through  the  elemental  boundaries  of  dr  and  is 
equal  to 


_0_  •  V^drdU 


II-9 


A  rigorous  derivation  of  Equation  II— 9  is  given  in  reference    [lj . 
Term  number  four  of  Equation  II— 4  is  the  absorption  term  and  accounts 
for  all  neutrons  lost  from  drdvd  XL     per  unit  time  by  absorption;    it 
is  given  by 


Z  ^  (  or)  (f>  Cr,  or,  xl, t )  d  r  d  ru-d _Tl 


11-10 
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The  rate  at  which  neutrons  are  scattered  out  of  drdvd-H-      into  new 
directions   jfl'     is  the  fifth  term  of  Equation  II-4  and  is  written  as 


Xs  CaO  ^Cr^XLjt  )  dr  dw'd-n-  .  11-11 


The  final  term  of  Equation  II— 4  is  the  time-rate  of  change  of  the 
neutron  density  in  drdvd-fl       .    This  term  can  be  written  as 


%i  n(r\  ru\  XL,t )  dr  d^dil  = 

TT      1  9 

if  §t  $cr,  n^-TL^t)  dr  d  r\r  d_n_   . 


By  placing  Equations  II— 5  through  11-12  into  the  verbal  Equation  II— 4 
and  dividing  out  the  differential  drdvdxL     the  result  is  the  Boltzmann 
transport  equation  for  the  distribution  of  neutrons  of  all  energies  in 
space  and  time 

"aF   It  ^Cr,/ir,ivb)  +  £tO)  0CF,n^,  JL,t) 
+  JrL-V  0Cf;ru-;_n_;t)  *  S  Cr,  at,  _fL}t)  11-13 

The  second  term  in  Equation  11-13  is  the  total  loss  of  neutrons  due  to 
scattering  out  to  new  spaces  and  directions  plus  absorption.     The  total 
cross  section        X-t  C  H>0        is  used. 

Equation  11-13  is  perfectly  general  with  the  following  two 
exceptions:    a)  the  medium  is  assumed  to  be  isotropic  and  homogeneous 
and  b)  neutron  slowing  down  is  due  only  to  elastic-scattering  collisions 
with  nuclei . 
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of  the  gth  group  is  defined  by  the  integral 

0gCF)  =  ffcr,  E)  dE  n-17 

where  the  upper  limit  on  the  integral    Eq_(       is  the  lower  boundary 
of  the  energy  interval  and  the  lower  limit  on  the  integral  Eg  is  the 
upper  boundary  of  the  energy  integral.    The  energy-dependent  flux  at 
the  point  F  is    ^CP^E") 

Within  each  group  the  diffusion  of  neutrons  is  described  by  the 
average  diffusion  coefficient.     For  simple  diffusion  theory  this  is 
written  as 


where    SI-(-r)q         is  the  group  macroscopic  transport  cross  section  and 
may  be  obtained  from 


where  the  summation  is  over  all  elemental  components  and  the  group 
microscopic  cross  section     <=r^.r  in  the  absence  of  pronounced 

resonance  effects  due  to  capture  and/or  scattering  with  the  energy 


group  is  given  by 


E3H 


z       _  J£SL    °"tZr)9CE)    0CE)d£ 


<=rtr     - ^r-  11-20 


'EgL 
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B.      THE  AGE-DIFFUSION  APPROXIMATION 

The  age-diffusion  approximation  to  Equation  11-13  is  developed 
in  reference  [l]  by  expanding  the  Boltzmann  transport  equation  in 
spherical  harmonics . 

In  the  steady  state  the  resulting  first  order  approximation  to 
Equation  11-13  is 


Since  the  change  in  the  slowing  down  density  of  neutrons  within  the  lethargy 
interval  du  must  be  balanced  by  the  flow  of  neutrons  scattered  into 
du  from  elsewhere,  the  last  term  of  Equation  11-14  may  be  written  as 


_B 


r 


And  Equation  11-14  may  be  rewritten  substituting  energy  for  lethargy 
variable  as 

11-16 

S<?,&)    -v-   1    Z5  CE.'-*  el")  0<O  d.E.' 

Complete  details  of  this  analysis  may  be  found  in  reference  [lj. 

C.      MULTIGROUP  FORMULATION 

In  forming  the  multigroup  diffusion  equation,  the  initial  step  is 
to  divide  the  complete  neutron  energy  (lethargy)  range  into  N  groups 
which  may  or  may  not  be  of  the  same  width.    The  neutron  flux 
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Neutrons  may  be  removed  from  a  particular  group  g  through  an 
absorption  interaction  described  by  the  group  absorption  cross  section 

^-  cxjQ   ~  —  r;Q  "**  — ^  c  >  q  +  ^-  m>0  "*"  -—  er-jQ        ;  II-21 


where 

Z.  r,q       =  fission  cross  section  of  group  g 
Xc  q        =  capture  cross  section  of  group  g 

^•m.3      =  Zv  "^ih,g^kc  =  inelastic  removal  cross  n_22 

^5*9  section  from  group  g 

^  enq     =  XT  ^  er.^K  =  elastic  removal  cross 
K^q  section  from  group  g 

The  group  macroscopic  cross  sections  defined  in  Equation  11-21  above 
are  obtained  from  the  corresponding  group  microscopic  cross  sections 
for  each  distinct  element  or  isotope  in  the  reactor  such  that 


,3  =  Z2   NZ   cr%)(5     7  11-23 


2 

where  q  is  the  particular  reaction  and  N     equals  the  number  of  atoms 
of  material  z  per  cubic  centimeter. 

The  multigroup  microscopic  cross  sections  for  each  material  are 
defined  in  reference  [7]  as  follows: 


GROUP  CAPTURE 
_  JE__ 11-24 

0(E)dE 

JEgL 
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GROUP  FISSION 


eqh 


cr 


5>9 


lL       J         r 

[^  $CE)dE 
Je9u 


11-25 


GROUP  SCATTERING  (IN  OR  OUT) 

CTZn,         (E-E')^fE)dEdE' 


m>^K 


rE^ 


toU 


?i(E)dE 


11-2  6 


where  E   TX  and  E   ,  are  the  upper  and  lower  limits  of  group  g  respectively 
gH  gL 

and     <fi  (E)  is  the  flux  per  unit  energy  interval  at  E. 

Neutrons  may  be  added  to  a  particular  group  g  by  a  fission  source 


N 


8 "' 


11-27 


where  %  g  is  the  fraction  of  fission  neutrons  born  into  group  g  such  that 

fE3H 
%3  =  %(E)dE         , 


11-28 


11-29 


and 


eooj)3  = 


rE3H 


£qL 


T)"CE)a-^f(E>)CZ)CE)dE 


^CE)dE 


11-30 
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The  steady  state  diffusion  equation  for  the  gth  energy  group,  where 
g  represents  any  group  between  one  and  N,  can  now  be  written 


D3V203(r)-Za,3  03(r) 
+    ^3  Z_.  ^h^f.h  ^h  (r)    =      S^Cr)         .  n-31 

In  Equation  11-31  the  inelastic  removal  cross  section  from  group  g 
appears  as  the  third  term  of  the  equation  rather  than  as  an  implied 
part  of  the  second  term.     This  equation  now  represents  a  set  of  N 
coupled  partial  differential  equations  independent  of  time  and  applicable 
to  one  region  of  a  reactor.    The  right  hand  side  of  the  equation 
represents  the  number  of  neutrons  at  energies  within  group  g  coming 
from  any  extraneous  neutron  sources. 

D.     SPACE-INDEPENDENT  FORM  OF  MULTIGROUP  EQUATION 

For  this  case,  it  is  assumed  that  the  extrapolation  distance  at  the 
surface  of  the  reactor  is  identical  for  all  energy  groups  and  thus  is 
independent  of  energy.     If  this  is  true  then  all  energy  groups  have  the 
same  spatial  dependence.     The  flux  is  then  written  as 

03(F)  '=    03  F(r)  11-32 

where     Gg  is  the  magnitude  of  the  flux  in  the  gth  energy  group,  and 
the  function  F(r)  satisfies  the  equation 

V2F(F)  +  BzF(r)  -O  11-33 

2 
the  factor  B     is  the  square  of  the  first  eigenvalue  of  Equation  11-33 

and  is  called  the  buckling. 
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By  substituting  Equations  11-32  and  11-33  into  Equation  11-31 
expanding  the  resulting  equation  in  a  series  of  eigenfunctions  ,  and 
applying  the  orthogonality  condition,  the  result  will  be  a  set  of  N 
linear  algebraic  equations  in  N  unknowns    $>  ?  $2)  '"   ^n       written 


9-1  h>3  N 


11-34 


This  set  of  equations  is  easily  solved  by  Cramer's  Rule. 

E.      EXPLICIT  SPATIAL  DEPENDENCE 

In  order  to  solve  the  set  of  Equations  11-31  by  numerical  methods, 
the  Laplacian  operator  in  the  first  term  must  be  replaced  by  its  equiva- 
lent finite  difference  approximation.     The  resulting  equation  for  a 
typical  energy  group  g  in  cylindrical  geometry  is 


11-35 


Appendix  A  contains  the  details  of  the  derivation  of  the  finite  difference 
approximation  and  the  definition  of  the  symbols  used  in  Equation  11-35. 
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Throughout  this  investigation,  circular  summetry  was  imposed. 
Thus  the  third  term  in  Equation  11-35  is  zero.     For  mesh  points  not 
on  the  longitudinal  axis  (See  Appendix  A  for  finite  difference  operator 
used  on  this  axis.),  Equation  11-35  becomes 

+  27h  I  ^9(r  +  l>z^  -  03<r-\,z:>3  +-RI  {  $g(r,  z +0 
-2  0g(r}z)  +  09Crlz-o]]  -2a,g  <^g  Cr,z.) 

*>S      N  h-, 


F.      MATRIX  FORM  OF  MULTIGROUP  EQUATIONS 

For  the  purpose  of  numerical  analysis  and  computational  ease, 
Equation  11-36  is  reduced  to  the  matrix  form.     In  Appendix  B,  this 
equation  is  first  expanded  to  fit  the  available  sixteen  group  data, 
then  reduced  to  a  form  more,  suitable  for  matrix  notation.     Finally, 
it  is  reduced  to  the  form 

A I  0g  (r-i,  z)+  A2  0g(r,z)  +  A30g  Cr+i.z^ 
+A409(r,z-D  +A5  0g(r) -z  +  0   =    Sg  (r)z) 

B-4 

+  AS,  0t  (r,z)-v  •  •   •    +AS9-1   (f)o)-\  (r,z) 
+  AF3  ^3  <r,z")  +  •  •  •  +  AFl8  0I6  C^z} 
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To  reduce  this  set  of  multigroup  equations  to  matrix  form,  the 
spatial  dependence  is  considered  first. 

The  flux  vector       <^>9         of  group  g  is  defined  as 


03    = 


11-37 


09<i>O 


0g(L,O 


9*9(1, M) 

09(L,m) 

where    ^gCL^^M*)      and     <^g  (2-»-L  3  M")        are  on  the  extrapolated 
boundary  and  are  set  equal  to  zero.     The  fluxes     0Ct>X)    and   0(x>\) 
are  set  equal  to  the  fluxes      0C3..X)      and        0("X>3)        respectively 
in  order  to  satisfy  the  normal  derivative  boundary  condition.     See 
Appendix  C   for  mesh  point  numbering  scheme. 
The  matrix  of  coefficients  Ag  is  written 

OAAO-'-O     AIA2A30--0     OA50--0--0 
O    O   A4  •  •  O     OAIA2A3--0     OOA50--0 


Aq   = 


o 


O      O 


A4      Q--AI 


A5 


•)  11-38 


The  terms  on  the  right  hand  side  of  equation  B-4  are  called  the 
source  terms  and  Wg  is  the  source  vector  for  group  g.    It  can  be 
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represented  in  matrix  form  by 


3-< 

z 


\Q 


W9    =  ]>T   ^Sh  0h    +  ^_,  AFh  0h  +    S9  0w 


where  the  flux  vector     9S  h      is  given  as 


<P 


K     = 


0hO,O 

0h(L,O 


0hCL;M) 


the  scattering  matrix  ASh 

O  •  •    •  O        O  ASk  O  •  -O         O 
O  •   -    -  O       OO     ASh  -Q        O 


ASu  = 


i_ 


o 


O  O 


the  fission  matrix  AFh 


AFh  - 


o-  - 

-    0 

OAFhO'    • 

O 

O 

0  ■  • 

•  0 

O    O  AFh- 

•O 

O 

o 


o      o 


o        o 


11-39 


11-40 


o 
o 


11-41 


O  O  •   •    '  ASy, 


O 
O 


11-42 


AF-, 
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and  the  extraneous  source  matrix  Sg 

o  -  •  •  -  o     o  s3  o  •  •  •  o    o 
0----0     oos^.-o     o 


o  • 


o    o 


o     o 


■o 

•o 


Sc 


11-43 


The  result  is  then  a  set  of  equations,  one  for  each  spatial  mesh 
point  for  the  energy  group  g.    These  may  be  written  as 


A3     03    =  Wg 


11-44 


Next  the  equations  are  cast  according  to  the  group  dependence. 
Define     c£>     as 


A  as 


4>  - 


L 


A,  O   • 
O    A2 


A  = 


0< 


o 
o 


11-45 


11-46 


o 


A< 
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and  W  as 


W 


w   = 


11-47 


wN 


The  entire  set  of  equations  then  takes  the  form 


A  0  =  W 


11-48 


G.      NUMERICAL  SOLUTION  OF  THE  MULTIGROUP  EQUATIONS 

The  basic  equation  11-48  can  be  solved  by  various  direct  or 
iterative  methods.    Direct  methods  generally  are  based  on  either 
Gaussian  elimination  with  pivoting,  or  on  triangular  decomposition  of 
the  matrix  of  coefficients.     Iterative  methods  include  the  Jacobi 
method,  unextrapolated  Liebmann  or  Gauss-Seidel  method,  various 
forms  of  the  extrapolated  Liebmann  or  successive  orerrelaxation 
method  and  others . 

In  this  study,  up  to  25,000  algebraic  equations  were  expected  as 
a  result  of  a  sixteen  group  problem  in  two  dimensions  with  1600  mesh 
points.     Tor  this  large  number  of  equations,  iterative  methods  are 
generally  used  rather  than  direct  methods  £8,  9].    Smith  [9]    states 
that  one  reason  for  this  is  that  the  iterative  methods  are  more  efficient 
in  that  they  take  advantage  of  the  sparse  matrix  of  coefficients  with 
its  numerous  zero  elements. 

In. this  study,  the  point  successive  overrelaxation  method  was 
used  (See  Appendix  D) . 
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First  the  reactor  materials  and  energy  group  structure  are  chosen. 
Next  the  cross  section  data  is  computed  or  as  was  the  case  in  this 
study,  taken  from  some  reliable  source    [7]  and  used  to  compute  the 
appropriate  coefficients.    Then  initial  estimates  of  the  group  fluxes 
designated       0g       at  all  mesh  points  in  the  reactor  except  those  lying 
on  the  extrapolated  boundaries  are  made.    These  values  are  then  used 

along  with  their  respective  coefficients  to  compute  the  initial  source 

—  ° 
vectors     Wo       .    The  flux  in  group  one  at  all  mesh  points  in  the 

reactor  is  then  corrected  by  solving  the  set  of  equations 


A,  0°    =  W,° 


11-49 


Using  the  point  successive  overrelaxation  method,  this  is  done 
by  solving  the  first  of  Equations  11-49  for  the  central  mesh  point, 
i.e.,      0,02,2)         and  then  solving  the  remaining  equations  point- by- 
point,  row-by-row,  at  all  other  points  no  on  the  extrapolated  boundary. 
The  solution  is  denoted    0,      .    The  remaining  group  equations  are  then 
similarly  solved  yielding  the  vector    0      where 


0'- 


0 


<Pu 


11-50 


which  represents  the  first  iteration. 

This  vector  is  then  rescaled  such  that 


<$,«*',  3.0')  =<  ¥°>  $°)   , 


11-51 
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which  indicates  that  the  flux  vector  on  the  zeroth  iteration  integrated 

over  all  spatial  dimensions  and  all  energy  groups  is  equal  to  the  flux 

vector  after  the  first  iteration  integrated  over  all  spatial  dimensions 

and  all  energy  groups  multiplied  by  a  scale  factor  [10]  .    When  the 

process  is  carried  out,  it  can  be  seen  that  after  a  few  iterations  the 

scale  factor  will  approach  a  constant  value.    The  reason  for  this 

as  described  by  Lamarsh  111]  is  that  each  iteration  is  actually 

equivalent  to  one  cycle  of  the  chain  reaction.    The  zeroth  iteration 

of  the  fluxes     <p      then  can  be  thought  of  as  a  set  of  initial  conditions 

or  the  flux  state  at  time  t  =  0.     Eventually  the  flux  will  approach  the 

funadmental  eigenfunction  as  the  higher  harmonics  die  out.    The 

fundamental  will  then  vary  with  the  number  of  iterations  (time) 

increasing  or' decreasing  depending  on  whether  the  system  is 

subcritical  or  critical.    The  quantity    ~£>      can  clearly  be  considered 

a  measure  of  the  criticality  of  the  system.     If     S     is  less  than  unity 

then  the  group  fluxes  throughout  the  system  are  decreasing  with  each 

iteration  and  the  system  is  subcritical  and  if    "5      is  greater  than 

unity  then  the  fluxes  are  increasing  and  the  system  is  supercritical. 

*£>      is  in  fact  the  inverse  of  the  effective  multiplication  factor  K   rr 

eff . 

The  iteration  process  described  above  is  referred  to  as  the  inner 
iteration  or  flux  iteration.    If  a  critical  assembly  is  desired,  then 
an  outer  iteration  must  be  performed  wherein  adjustments  are  made  to 
the  composition  or  size  of  the  reactor  in  order  to  make  the  effective 
multiplication  factor  converge  to  unity. 

Convergence  of  the  flux  vectors  on  the  inner  or  flux  iterations 
is  determined  by  taking  the  absolute  value  of  the  difference  between 
the  flux  at  all  mesh  points  and  in  all  energy  groups  and  comparing  it 
with  some  desired  value      £  (     ,  such  that 


<   €,  11-52 


I 

ALL 
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Convergence  of  the  outer  or  criticality  iteration  is  determined  by- 
subtracting  the  scale  factor    t>      from  unity  so  that 


i.o  -  s 


<  e 


z 


11-53 


where    £z  is  an  arbitrary  convergence  criterion. 
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III.     DESCRIPTION  OF  COMPUTER  PROGRAM 

A.  INPUT  DATA 

In  this  section  the  input  data  is  read  in  via  punched  cards  as 
illustrated  in  Appendix  E  and  printed  in  the  output  for  verification. 
The  input  data  along  with  their  symbolic  program  designators  include 
the  following:    the  energy  dependent  source  SS  (I);  the  maximum  values 
of  the  initial  flux  estimate  S(I);  the  energy  intervals  DELE(I);  the  number 
of  atoms  per  cubic  centimeter  H(I)  of  each  element  or  isotope;  the 
microscopic  cross  sections  for  each  element  or  isotope  for  inelastic 
scattering  SIGIN(I,J,K) ,  transport  SIGTR(IJ),  fission  SIGF(I  J), 
elastic  removal  SIGER(I,J),  and  capture  SIGC(IJ);  the  normalized 
neutron  fission  spectrum  CHI(I);  and  the  average  number  of  neutrons  per 
fission  NU(IJ). 

B.  OUTER  ITERATION 

If  a  critical  composition  is  the  desired  result,  the  program  returns 
to  this  section  after  the  iteration  has  converged  (See  SAMPLE  COMPUTER 
PROGRAM  page  119).    The  iteration  on  composition  then  begins  by 
adding  to  the  volume  of  the  reactor  an  incremental  volume  of  uranium 
235.    An  equal  volume  of  uranium  238  is  subtracted  and  new  values  for 
N2^   H(l)     and  N28   H(2)     are  computed.     On  the  zeroth  outer  iteration 
this  section  of  the  program  is  by-passed. 

C.  COMPUTATION  OF  MATRIX  OF  COEFFICIENTS 

Input  data  is  used  in  this  section  to  compute  the  elements  of  the 
coefficient  matrix,  the  scattering  matrix,  and  the  fission  source  matrix. 
In  the  first  part  of  this  section,  the  various  microscopic  cross  sections 


for  each  reactor  component  are  multiplied  by  their  respective  atomic 
abundances  N    .    The  values  for  each  particular  reaction  (for  example 
capture  or  inelastic  scatter  from  g  to  g  +  1)  are  then  summed.     This 
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is  done  for  each  energy  group.     The  summations  are  then  combined 
with  other  input  data  and  the  coefficients  are  computed  for  each 
group  equation. 

The  equation  denoted  by  statement  number  17  in  the  sample 
program  computes  the  diffusion  coefficients  for  each  group  based  on 
the  simple  diffusion  theory  approximation.    The  absorption  coefficients 
for  each  group  are  computed  in  the  equation  denoted  by  statement 
number  175  by  summing  the  contributions  from  all  the  absorption 
or  loss  reactions.    The  equation  denoted  by  statement  number  18 
computes  the  fission  and  down-scattering  coefficients  from  all 
higher  energy  groups .     This  is  accomplished  by  starting  at  the 
highest  energy  group  in  which  the  reactions  originate  and  computing 
its  contribution  to  all  lower  energy  groups  at  the  end  of  the  slowing- 
down  process.    The  fission  source  coefficients  for  a  particular  group 
and  all  lower  energy  groups  are  computed  by  the  equation  given  in 
statement  number  21.    The  above  coefficients  are  completely  defined  in 
Section  II  and  Appendix  B. 

D.      INITIAL  FLUX  ESTIMATE  AND  BOUNDARY  CONDITIONS 

In  the  cylindrical  geometry  the  flux  shape  is  initially  estimated 
as  the  product  of  the  relative  group  magnitude  factor  S(I),  a  zero 
order  Bessel  function,  and  a  cosine  function.    The  factors  S(I)  were 
determined  for  this  study  as  suggested  by  Clark  and  Hansen   [10]  . 

The  flux  on  the  exterior  boundaries  is  set  equal  to  zero.    At  the 
interior  boundaries  the  flux  at  the  first  mesh  points  on  either  side  of 
the  longitudinal  and  mid-plane  axes  are  set  equal,  indicating  that  the 
flux  has  its  maximum  value  along  these  axes.    At  the  end  of  this  section 
the  initial  flux  estimate  is  printed  out  for  verification  and  comparison 
with  the  final  flux  values . 

The  program  is  designed  to  by-pass  this  section  on  all  outer 
iterations . 
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E.         INITIAL  FISSION  SOURCE  INTEGRATION 

The  initial  flux  estimate  is  used  in  this  section  to  compute  the 
fission  source  vector  for  the  zeroth  inner  iteration.     The  program 
performs  the  following  triple  integration: 


z  re 


o  0 


{[^(E;r,z)]  X(E)^/l)>cCE)Zf,jL(E)jdEdrdz         m-i 


by  dividing  the  energy  spectrum  into  discrete  intervals  DELE  (I) 
corresponding  to  the  energy  groups  (for  this  study  the  upper  boundary 
of  the  highest  energy  level  was  taken  as  10  Mev  for  use  in  this  calcu- 
lation) and  dividing  the  area  into  finite  intervals  corresponding  to  the 
mesh  spacing  GS .    A  summation  process  is  then  carried  out. 

F.  COMPUTATION  OF  SOURCE  VECTORS 

The  components  of  the  vector  W  for  each  inner  iteration  are  computed 
by  summing  the  products  of  the  fission  source  or  the  combined  fission 
source  plus  down-scattering  source  coefficients  and  the  scaled  flux 
values  from  the  previous  iteration.     On  the  zeroth  iteration  the  initial 
flux  estimate  is  used  with  a  scale  factor  of  unity.     The  resulting  numbers 
SUMR(I,R,Z)  and  SUML(I,R,Z)  are  the  space  dependent  components  of 
each  group  source  vector  Wg. 

G.  SOLUTION  OF  THE  DIFFUSION  EQUATION 

This  section  consists  of  one  outer  energy  group  loop  with  several 
space  loops  within  it.     Each  group  equation  beginning  with  the  highest 
group  is  solved  at  all  spatial  points  starting  at  the  longitudinal  axis  r  =  0 
and  proceeding  column  by  column  to  within  one  mesh  point  of  the  extra- 
polated boundaries  r  =  R-l  and  z  =  Z- 1 . 
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The  first  space  loop  will  solve  the  diffusion  equation  at  the 
geometric  center  of  the  reactor,  i.e.  ,  point  (2,2)  (See  Appendix  C, 
Figure  C-l.)  if  an  extraneous  point  source  is  located  in  this  position. 
This  particular  configuration  was  used  in  this  study;  however,  extrane- 
ous sources  with  other  spatial  distributions  could  also  be  used  with 
the  appropriate  changes  to  this  section  of  the  program. 

The  second  space  loop  involves  the  solution  of  the  diffusion 
equation  at  all  points  on  the  longitudinal  axis  r  =  0  but  not  occupied 
by  an  extraneous  source.    A  special  finite  difference  approximation 
for  points  on  this  axis  is  necessary  due  to  a  singularity  caused  by  the 
radius  being  equal  to  zero  (See  Appendix  A,  page  66). 

The  last  spatial  loop  solves  the  diffusion  equation  at  all  remaining 
mesh  points  with  the  exception  of  those  located  on  the  exterior 
boundaries  . 

Within  each  of  these  spatial  loops,  the  flux  value  at  each  mesh 
point  from  the  previous  iteration  is  unsealed  (TPHI1)  and  compared 
with  the  flux  value  at  each  mesh  point  in  the  present  iteration  to 
determine  convergence  as  defined  in  Section  II  ,    page  27.    The 
largest  absolute  value  of  these  differences, 

TEST2     =    DABS  (TPH1Z  -  TRUll}     ,  III-2 

is  stored  as  the  quantity  SAVE  and  may  be  printed  with  the  output  as 
desired.    The  scaled  flux  values  of  the  previous  iteration  SPHI1, 
SPHI2,  SPHI3  are  used  in  the  successive  point  overrelaxation  scheme 
for  solving  the  group  equations  at  each  mesh  point. 

H.      NTH  FISSION  SOURCE  INTEGRATION  (UNSCALED  FLUX) 

In  this  section  the  triple  integration  process,  Equation  III—  1 ,  is 
performed  with  the  new  but  as  yet  unsealed  flux  values. 
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I.       COMPUTE  SCALE  FACTOR  AND  RENORMALTZE  FLUX 

Once  the  fission  source  integration  has  been  computed  for  the 
N-l  iteration  FLUXI(NN-l)  using  the  scaled  or  renormalized  flux 
values  and  for  the  Nth  iteration  FLUXI(NN)  using  unsealed  flux  values, 
the  scale  may  be  determined  from 

SCALF  =  DSQRT(FLUXI(NN-1)/FLUXI(NN))  III- 2 

The  flux  values  at  all  mesh  points  and  in  all  energy  groups  are  then 
renormalized  by  multiplying  by  this  scale  factor.    The  scale  factor 
SCALF  is  the  principal  eigenvalue  of  the  system  and  its  reciprocal  is 
the  effective  multiplication  constant  SKEFF  as  explained  in  Section  II, 
page 

J.       NTH  FISSION  SOURCE  INTEGRATION  (SCALED  FLUX) 

At  this  point  the  integration  process  of  Section  III (E)  is  performed 
using  the  renormalized  flux.    On  the  subsequent  iteration  this  will 
become  the  quantity  FLUXI(NN-l)  in  the  numerator  of  Equation  III— 2 . 

K.      END  OF  INNER  ITERATION 

This  statement  (1801)  is  the  controlling  statement  for  the  inner 
iteration.     It  can  make  the  program  exit  to  either  the  output  section 
or  the  outer  iteration  loop  depending  upon  whether  or  not  a  critical 
reactor  is  desired.    The  controlling  factor  may  be  either  a  specific 
number  of  inner  iterations  if  experience  has  shown  that  the  convergence 
criterion  as  defined  in  Section  II,  page  27  will  be  satisfied  or  another 
possibility  would  be  to  stop  the  inner  iteration  once  the  convergence 
criterion  is  met  regardless  of  the  number  of  inner  iterations.    The 
second  method  might  be  preferable  once  it  had  been  determined  that 
convergence  to  a  reasonable  criterion,  i.e.  ,  6,  would  not  take  an 
excessive  amount  of  computer  time. 
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L.      OUTPUT  SECTION 

The  output  is  very  flexible  and  may  be  easily  changed  to  fit  a 
particular  problem  or  set  of  problems.    Typical  output  data  includes  the 
following:    scale  factor,  effective  multiplication  constant,  flux 
values  at  various  mesh  points  and  various  inner  or  outer  iterations, 
convergence  criterion  SAVE  (  6j  )  and  overrelaxation  factor  w.    In 
addition,  graphical  output  in  the  form  of  graphs  of  flux  shape  compared 
with  the  cosine  function  and  zero  order  Bessel  function  at  various 
locations  may  be  plotted  and  the  convergence  criterion  SAVE  versus 
iterations.     Many  other  possibilities  exist. 

M .     START  NTH  OUTER  ITERATION  AND  TEST  FOR  CRITICALITY 

In  this  section  the  scale  factor  is  checked  out  for  criticality 
by  comparison  with  an  arbitrarily  selected  value  (  .  0005  was  chosen 
for  this  study.)  such  that 

|  1.0  -  SCALF  1        £      .0005  III-3 

If  the  above  condition  is  not  satisfied,  then  another  outer  iteration 
is  made.     The  number  of  outer  iterations  may  also  be  arbitrarily 
limited  in  this  section.    Once  this  limit,  for  example,  ten  outer 
iterations,  is  reached,  the  program  will  stop  even  though  criticality 
has  not  been  reached. 
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IV.     RESULTS  AND  PROGRAM  VERIFICATION 

If  in  the  paragraphs  that  follow  modifications  to  the  basic  core 
described  below  have  been  made,  these  modifications  will  be  so  stated 
The  results  described  and  tabulated  in  this  section  were  based  on  a 
typical  homogeneous  fast  reactor  core  composed  of  the  following 
materials  with  their  respective  volume  fractions: 

U235-      2"/ 
U238  -  30% 

Na       -  33% 

Fe239"  27% 
Pu        -     8% 

This  composition  was  chosen  to  coincide  approximately  with  the 
RAPSODIE  reactor  core  in  reference  [  12]  .     The  core  dimensions  were 
36  centimeters  in  radius  and  72  centimeters  in  height.     Mesh  spacing 
was  two  centimeters  in  both  the  radial  and  axial  directions.    Three 
energy  groups  were  used  with  the  cross  section  data  taken  from 
ANL  5800    [7]. 

A.      DETERMINATION  OF  OPTIMUM  RELAXATION  FACTOR 

The  experimental  method  used  for  determining  w         consisted  of 

opt 

two  basic  steps.     First,  an- analytical  estimation  was  obtained  based 
on  formulation  developed  by  Frankel  [13]  and  Starling  £14]  . 

Frankel   £13]  derived  the  following  equations  based  on  the 
solution  to  Laplace's  equation  in  a  bounded  region  R  in  rectangular 
coordinates 


wopt 


=   A  OC  IV- 1 
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where   ©c  is  the  smaller  root  of 

Oc2t2   -  A  ol    -t-   I      -O  ,  IV-2 

and  t  is  given  by 

t    =     cos^-     +    COS    X  IV- 3 

where  p  and  q  are  the  number  of  mesh  points  in  the  x  and  y  directions. 
Smith  £9j  used  this  same  method  to  solve  the  similar  Poisson's  equa- 
tion in  rectangular  coordinates  in  a  bounded  region  R.     Starling  £14] 
has  developed  an  equation  for  t  applicable  to  Laplace's  equation  in 
cylindrical  coordinates  and  given  by 

t    =    cos  ^     +   To  (       R^.—  )        ,  rv-4 

where  Z  and  R  are  the  height  and  radius  of  the  oyJinder  respectively 

and  h  is  the  mesh  spacing. 

The  table  below  gives  the  results  found  for  these  analytical 

solutions  for  w        using  a  symmetric  cylindrical  section  36  centimeters 
opt 

in  radius  and  36  centimeters  in  half-height 

METHOD  FOR  FINDING  "t"  w      t  REMARKS 

opt 

Eq.  IV- 4  1.7888  Z  =  72  cm .   (full  height) 

Eq.  IV- 4  1.7259  Z  =  36  cm.   (half  height) 

Eq.  IV- 3  1.7294  p  =  q  =  20  (no.  mesh  pts .  +  1) 

Eq.  IV- 3  1.704  p  =  q  =  18  (no.  mesh  pts.) 

The  second  step  in  the  experimental  process  was  to  try  these  values 

in  the  actual  program  and  graph  the  convergence  criterion  versus  the 

number  of  inner  iterations  for  these  values  of  w         along  with  other 

opt 

reasonable  estimates  to  determine  the  best  one.     Figure  IV- 1  shows  the 
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data  which  resulted  from  this  study.     It  can  be  seen  that  for  the  values 

of  w         equal  to  1.72,    1.73  and  1.74,  a  convergence  criterion  value 

op  -6 

of  less  than  10       is  satisfied  within  fifty  iterations.     Based  on  this 

data,  it  was  concluded  that  Equation  IV- 4  with  Z  equal  to  one-half 

the  height  of  the  cylinder,  due  to  symmetry,  gave  the  best  value  of 

w        for  this  program.    The  program  was  modified  to  compute  w 
opt  opt 

automatically  based  on  Equation  IV- 4. 

The  above  results  were  reasonable  in  view  of  the  fact  that 
Equation  IV- 4  was  developed  for  use  with  Laplace's  equation  in 
cylindrical  coordinates .    The  multigroup  neutron  diffusion  equation 
used  in  this  study  in  cylindrical  geometry  is  of  a  similar  form.     One 
of  the  major  differences  is  that  the  dependent  variable   0     is  a 
function  of  both  space  and  energy . 

B.      FLUX  SHAPES  FOR  A  TYPICAL  NONCRITICAL  CORE 

For  the  bare  reactor  described  at  the  beginning  of  this  section,  the 
initial  flux  shape  was  estimated  to  be  the  product  of  the  relative  group 
weighting  factor  S(I),  a  cosine  function  and  a  zero  order  ordinary 
Bessel  function  of  the  first  kind.     This  corresponds  to  the  solution  of 
Laplace's  equation  in  a  finite  cylinder.     Figure  IV- 2  shows  the  final 
flux  shape  in  the  axial  direction  at  the  center  of  the  reactor  after 
260  inner  iterations  with  a  convergence  criterion  of  less  than  10 
compared  with  the  cosine  shape  at  the  same  location  at  zero  iterations. 
The  magnitudes  were  normalized  to  unity  so  that  the  shapes  could  be 
easily  compared.    It  was  found  that  the  deviations  were  in  the  order  of 
1  to  15  per  cent  with  the  per  cent  deviation  increasing  from  the  center 
outward.    All  three  energy  groups  showed  exactly  the  same  normalized 
shapes  after  an  equal  number  of  iterations.    The  magnitudes  of  the 
group  fluxes  are  unknown  since  they  depend  on  the  power  level  at 
which  the  reactor  is  being  operated. 


39 


AXIAL  FLUX  SHAPE 
IN  A  TYPICAL  SUBCRITICAL  CORE 
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The  same  comparison  has  been  made  in  Figure  IV- 3  between  the 
flux  in  the  radial  direction  and  the  assumed  initial  Bessel  function. 
A  very  slight  deviation  in  shapes  was  noted  here  also  indicating  the 
close  similarity  between  the  diffusion  equation  and  Laplace's 
equation. 

Examination  of  the  data  revealed  that  the  final  flux  shapes  were 
in  effect  established  completely  to  seven  decimal  places  after  some  60 

inner  iterations  at  which  time  the  convergence  criterion  SAVE  was  less 

-7 
than  10      . 

The  program  was  run  using  a  constant  as  the  initial  flux  shape 

estimate  rather  than  the  cosine  and  Bessel  function.     The  resulting 

final  shapes  for  all  three  groups  were  exactly  the  same  as  those  in 

Figures  IV -2  and  TV- 3  after  90  inner  iterations.    Convergence  was 

somewhat  slower  using  this  initial  estimate.     For  example,  after  50 

-5  -7 

iterations,  SAVE  was  equal  to  3  . 2  X  10       versus  3.7  X  10       for  an 

initial  estimate  based  on  the  cosine,  Bessel  function  product. 

C.      PROGRAM  TEST  WITH  EXTRANEOUS  POINT  SOURCE 

The  program  was  run  with  a  simulated  extraneous  Pu-Be  point 
neutron  source  at  the  geometric  center  of  the  reactor.    Two  core 
compositions  were  used  in  these  tests. 

The  first  core  was  made  up  of  the  typical  composition  described 
at  the  beginning  of  this  section.     The  flux  shapes  in  the  radial  and 
axial  directions  for  the  three  groups  were  plotted  on  the  same  graph 
with  the  Bessel  functions  and  cosine  respectively.     Figures  IV- 4  and 
rV-5  show  the  results  of  these  comparisons  after  90  inner  iterations. 
Group  one  showed  the  largest  deviation  with  the  most  pronounced 
effect  nearest  the  source.    This  result  was  due  to  the  normalized  source 
neutron  contribution  to  group  one.    Eighty-seven  per  cent  of  the  neutrons 
from  the  Pu-Be  source  are  born  with  energies  lying  in  group  one 
(1.35  Mev  ->  co    )  as  defined  in  ANL  5800  [  7]. 
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RADIAL  FLUX  SHAPE 
IN  A  TYPICAL  SUBCRITICAL  CORE 
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RADIAL  FLUX  SHAPES 

IN  A  TYPICAL  SUBCRITICAL  CORE 

WITH  A  CENTRAL  POINT  SOURCE 
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AXIAL  FLUX  SHAPES 
IN  A  TYPICAL  SUBCRITICAL  CORE 
WITH  A  CENTRAL  POINT  SOURCE 
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The  second  core  was  a  solid  iron  cylinder  with  the  same 
dimensions  as  the  typical  fast  reactor  core.     The  Pu-Be  source 
was  located  at  the  center  of  the  cylinder.     The  radial  flux  shapes  for 
groups  one  and  two  were  plotted  in  Figure  IV-6  and  the  axial  flux 
shapes  for  these  two  groups  were  plotted  in  Figure  IV- 7.    It  was  noted 
that  for  these  two  groups  the  flux  decayed  very  rapidly  as  the  distance 
from  the  source  increased;  therefore,  the  ordinates  of  Figures  IV-6 
and  IV- 7  are  logarithmic.    Figures  IV- 8  and  IV-9  show  the  radial  and 
axial  flux  for  group  three.    In  these  graphs  the  flux  shapes  are  com- 
pared with  the  cosine  and  Bessel  functions  after  90  inner  iterations. 

The  three  group  flux  spectrum  at  the  center  of  the  two  assemblies 
are  shown  in  Figure  IV- 10.    A  relatively  flat  spectrum  was  observed  for 
the  core  containing  the  fissile  material  while  the  solid  iron  core  showed 
large  relative  differences  among  the  three  group  fluxes. 

u  . 

Test  of  the  program  up  to  this  point  included  only  the  inner  iteration. 
Criticality  was  achieved  by  two  methods  of  outer  iteration. 

The  first  method  consisted  of  performing  sufficient  inner  iterations 

for  convergence.    Seventy  iterations  were  used  resulting  in  an    £t     of 

-9 
less  than  10       for  a  mesh  size  of  two  centimeters.    An  outer  iteration 

235 
was  then  carried  out  wherein  .5%  by  volume  of  U         was  added  before 

each  outer  iteration  as  described  in  Section  III-B,    This  cycle  was  then 

repeated  until  the  criticality  criterion    (       |  1.0  -  5KEFF  \<.Ooo5  ) 

was  satisfied.    The  optimum  overrelaxation  factor  remains  the  same 

for  each  outer  iteration. 

Four  different  mesh  spacings  were  used  in  this  section  of  the 

program  analysis  to  determine  the  effect,  if  any,  mesh  size  had  on 

criticality  and  flux  shapes.    Table  IV- 1  contains  the  data  resulting 

from  this  study.    The  central  flux  shapes  in  the  radial  and  axial 

directions  at  criticality  for  all  mesh  sizes  tried  were  nearly  coincident 
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AXIAL  FLUX  SHAPE 

IN  AN  IRON  CYLINDER 

WITH  A  CENTRAL  POINT  SOURCE 
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RADIAL  FLUX  SHAPE 

IN  AN  IRON  CYLINDER 
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FLUX  SPECTRUM 
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WITH  A  CENTRAL  POINT  SOURCE 
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with  the  Bessel  function  and  cosine  respectively.     Table  IV- 2  shows  a 

maximum  deviation  of  .71  per  cent  using  a  mesh  size  of  two  centimeters. 

235 
Criticality  was  reached  at  the  same  core  composition,  i.e.  ,  4.5%  U        , 

for  mesh  spacings  of  .75  centimeters,    1.5  centimeters  and  2  centimeters. 

With  a  mesh  size  of  four  centimeters,  criticality  was  not  achieved 

at  this  composition. 

The  second  method  for  arriving  at  a  critical  assembly  consisted 

of  fixing  the  core  composition  and  varying  the  core  dimensions  in  the 

outer  iteration.    A  new  optimum  overrelaxation  factor  was  computed 

for  each  outer  iteration  since  the  geometric  parameters  that  were  used 

to  compute  w        were  changing.    Three  compositions  corresponding 
opt 

to  18.9%,   20%  and  27%  enrichment  were  used.     The  results  of  this 
test  are  plotted  in  Figure  IV- 11  and  compared  with  the  bare  core  criti- 
cality curve  of  the  RAPSODIE  reactor  in  reference  [11]  .    Differences 
noted  are  believed  to  be  due  to  the  following:    a)  compositions  are  not 
exactly  the  same;  b)  different  cross  section  data  was  probably  used 
though  this  was  speculative  since  reference  [11]  did  not  state  the 
source  of  their  cross  section  data;  and  c)  different  programs  were  used 
to  solve  the  multigroup  diffusion  equations.     PRODII  program  was  used 
for  the  RAPSODIE  reactor.    The  factor  above  having  the  most  influence 
was  probably  the  large  volume  of  iron.     The  absorption  cross  section 
of  iron  is  fairly  large  and  as  seen  in  Table  IV- 3  is  occupying  a  volume 
which  in  RAPSODIE  contained  4%  molybdenum,  an  unspecified  quantity 
of  iron  in  the  form  of  stainless  steel,  and  possibly  some  void  space. 
The  absorption  cross  section  of  these  components  is  likely  to  be  less 
than  pure  iron . 

The  comparison  of  these  curves  shows  that    while  a  much  larger 
critical  reactor, assuming  the  same  enrichment,  would  result  from 
using  the  program  developed  in  this  study,  the  shapes  of  the  curves 
are  very  similar.    The  fact  that  as  enrichment  was  increased,  the  core 
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RAPSODIE 
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3 

3 

CROSS  SECTION 
DATA 

? 

ANL  3  GROUP 

ANL  3  GROUP 

COMPOSITION 
VOL.   PERCENT 

URANIUM 

28% 

32% 

32% 

PLUTONIUM 

8% 

8% 

8% 

SODIUM 

33% 

33% 

33% 

MOLYBDENUM 

4% 





STAINLESS  STEEL 
IRON 

? 

27% 

27% 

Table  IV- 3 
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size  and  hence  the  critical  mass  was  decreased  indicates  a  correct 
trend.    The  flux  shapes  found  using  this  outer  iteration  scheme  were 
also  nearly  coincident  with  the  cosine  and  Bessel  function  curves. 
Table  IV- 4  shows  an  example  of  the  per  cent  deviations  experienced 
for  a  critical  core  assembly  32  centimeters  in  radius  and  124  centimeters 
in  height. 

E.         TESTS  WITH  SIXTEEN  GROUP  DATA 

Attempts  to  modify  the  program  to  solve  a  sixteen  group  problem 
were  not  successful.     Both  under-  and  over-relaxation  factors  were 
tried  but  satisfactory  convergence  of  the  inner  iteration  could  not 
be  achieved  after  over  thirty  minutes  of  computer  time.    An  under- 
relaxation  factor  of  .  1  was  used  and  after  310  iterations  and  34 
minutes,  the  convergence  criterion  was  only  of  the  order  10 
Overrelaxing  with  a  factor  of  1.726  yielded  a  highly  oscillatory  flux 
behavior  in  all  energy  groups.    All  group  fluxes  were  initially  estimated 
to  have  maximum  value  of  ten.    The  magnitude  of  oscillations  recorded 
in  the  output  varied  between  plus  or  minus  eighty  for  the  higher  energy 
groups  to  plus  or  minus  thirty  for  the  lower  energy  groups  with  no 
noticeable  damping  after  310  iterations  and  some  34  minutes  of 
computer  time.     No  further  tests  were  conducted. 
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MESH  POINT 

INITIAL  FLUX 

FINAL  FLUX 

PER  CENT 

LOCATION 

GROUPS  1,2,3 

GROUPS  1,2,3 

DEVIATION 

(2,2) 

1.000000 

1.000000 



(4,2) 

.977535 

.977479 

.006 

(6,2) 

.911646 

.911554 

.010 

(8,2) 

.806757 

.806654 

.013 

(10,2) 

.669888 

.669786 

.015 

(12,2) 

.510080 

.510013 

.013 

(14,2) 

.337808 

.337782 

.008 

(16,2) 

.164121 

. 164154 

.020 

(18,2) 

0.0 

0.0 

• 

In     n\ 

1.000000 

1. 000000 



(2,4) 

.994868 

.994805 

.006 

(2,6) 

.979525 

.979403 

.012 

(2,8) 

.954141 

.953952 

.020 

(2,10) 

.918959 

.918717 

.026 

(2,12) 

.874343 

.874059 

.032 

(2,14) 

.820767 

.820438 

.040 

(2,16) 

.758757 

.758407 

.046 

(2,18) 

.688969 

.688603 

.053 

(2,20) 

.612101 

.611742 

.058 

(2,22) 

.528959 

.528615 

.065 

(2,24) 

.440393 

.440076 

.072 

(2,26) 

.347303 

.347033 

.078 

(2,28) 

.250656 

.250440 

.086 

(2,30) 

.151424 

.151289 

.089 

(2,32) 

.050646 

.050500 

.091 

Table  IV-4 
57 


V.     CONCLUSIONS 

This  study  was  undertaken  to  develop  an  easy  to  apply  computer 
code  for  numerically  solving  the  multigroup  diffusion  equations.     The 
program  was  written  in  the  FORTRAN  IV  language  and  all  work  was  done 
on  an  IBM  360/67  computer.    Applicability  to  other  computer  systems, 
while  not  a  consideration  of  this  study,  is  assumed  to  be  practical. 

The  unique  feature  of  this  program  is  that  it  will  be  readily 
available  to  and  easily  used  by  students  at  the  Naval  Postgraduate 
School.     It  is  hoped  that  this  study  will  complement  future  local 
experimental  studies  in  fast  neutron  diffusion  and  also  enable  the 
student  researchers  to  conserve  valuable  time  which  would  be 
expended  attempting  to  verify  their  experimental  data  with  a  myriad 
of  complicated  and  unfamiliar  multigroup  diffusion  programs. 

The  point  successive  overrelaxation  method  of  solution  was 
successful  with  the  three  group  da  La  on  the  basis  of  the  convergence 
criterion  chosen,  i.e., 

|0N-'CE)r>z)~  0NCE,r,z)l      <   \0'& 

ALL 

Convergence  of  the  three  group  criticality  problem  with  five  outer 
iterations  took  five  minutes  and  forty  seconds  and  126K- words  of 
core  storage  using  a  mesh  size  of  two  centimeters.     Criticality  was 
achieved  with  smaller  mesh  sizes,  i.e.  ,    .75  centimeters  and  1.5 
centimeters,  but  computation  times  were  considered  excessive  especially 
for  .75  centimeter  mesh  size  which  took  about  33  minutes  of  computer 
time  and  over  twice  as  much  core  storage  space. 

Tests  showed  that  the  code  was  easily  adaptable  to  three  group 
solutions  for  typical  fast  reactor  cores  and  reactor  material  assemblies 
with  extraneous  point  sources.     Flux  shapes  and  energy  spectrums  from 
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these  tests  should  provide  a  sound  basis  for  future  experimental  flux 
measuring  studies  in  fast  reactor  materials. 

The  close  coincidence  of  the  critical  flux  shapes  with  the  cosine 
and  Bessel  functions  indicates  that  for  a  three  group  problem  if  critical 
core  size  and  composition  are  known,  then  the  space-independent  solu- 
tion using  the  buckling  may  be  used.     It  was  found,  however,  that 
this  was  not  true  for  subcritical  assemblies,  i.e.  ,  the  space  and 
energy  dependence  of  the  flux  are  not  separable. 

Further  work  with  larger  numbers  of  energy  groups  will  have  to  be 
carried  out  to  establish  the  program's  ability  to  handle  successfully  many 
group  calculations.    The  unsuccessful  attempt  at  using  sixteen  groups 
was  due  to  the  failure  to  attain  suitable  convergence  in  a  reasonable 
amount  of  time.    A  better  technique  than  the  point  successive  over- 
relaxation  method  is  obviously  needed  when  dealing  with  this  number 
of  energy  groups  in  order  to  improve  the  convergence  time. 
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VI.     RECOMMENDATIONS  FOR  FUTURE  WORK 

In  its  present  form  the  program  has  proven  to  be  successful  for 
the  purpose  for  which  it  was  intended.    Future  researchers  might  choose 
to  expand  or  modify  this  program  in  one  or  more  of  the  following  areas. 

A.  GEOMETRY  AND  SPATIAL  DIMENSIONS 

Although  the  program  is  written  for  the  cylindrical  geometry  in  two 
spatial  dimensions  (flux     d      is  not  a  function  of  azimuth  angle    ©       ), 
it  need  not  be  so  limited.    The  program  could  easily  be  adapted  to  the 
three  dimensional  cylindrical  geometry  as  well  as  the  two  or  three 
dimensional  rectangular  geometry.    The  addition  of  a  third  dimension  in 
space  would,  of  course,  require  considerably  more  computer  core  storage 
and  hence  the  future  user  might  have  to  limit  the  physical  size  and 
number  of  energy  groups  in  a  particular  problem. 

Another  possible  additional  feature  might  be  to  increase  the 
number  of  regions  by  one  or  more  to  include  blanket,  reflector  and 
shielding  regions. 

B.  CROSS  SECTION  DATA 

All  the  results  for  this  program  were  obtained  using  three  group 
ANL  data.    Conceivably  future  users  would  want  to  use  the  most 
current  and  widely  accepted  cross  section  available  to  them.     For 
this  reason  the  program  could  be  modified  to  handle  multigroup  data  of 
up  to  twenty-six  groups  or  higher.     Once  again  the  user  must  bear  in 
mind  that  as  the  number  of  energy  groups  is  increased,  core  storage  and 
computation  time  will  also  increase  and  probably  not  linearly. 

C  .         INNER  ITERATION  METHODS 

The  point  successive  overrelaxation  iterative  method  gave 
satisfactory  convergence  rates  for  the  three  group  problems  used  in 
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testing  the  program.     However,  more  rapid  convergence  of  the  inner 
iteration  is  possible  [15]    using  line  rather  than  point  successive 
overrelaxation  methods  and  cyclic  Chebyshev  semi-iterative  methods. 
If  the  future  user  desires  to  use  more  than  three  energy  groups,  he  will 
almost  certainly  have  to  change  the  program  to  fit  one  of  these  methods 
of  solution, 

D.         TIME  DEPENDENCE 

On  a  somewhat  more  far  reaching  level,  this  program  might  be 
expanded  to  include  the  time  domain  for  the  purpose  of  studying  the 
dynamics  of  a  fast  reactor.    This  could  also  include  a  study  of  space- 
dependent  burnup  and  reactivity  feedback  effects.     These  topics  are 
currently  under  investigation  by  those  working  on  the  development  of 
the  fast  breeder  reactor. 
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APPENDIX  A 

DERIVATION  OF  LAPIACIAN  OPERATOR  IN  CYLINDRICAL 

COORDINATES  AND  FINITE  DIFFERENCE  APPROXIMATIONS 


In  order  to  numerically  solve  the  set  of  partial  differential 
equations  11-31,  a  finite  difference  technique  was  chosen.     In 
cylindrical  geometry  the  Laplacian  or     V        operator  is  approximated 
by  transforming  the  Cartesian  coordinates  to  polar  coordinates. 


X 


P(x,y,z.)  =  PCr,^,©) 


x  -  r cos e 


y  =    r  sine 


r  =  cx2+y2)/2-    ■      Gstah-'^ 


A-l 


> 
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Taking  the  partial  derivations  of  r  and   0     with  respect  to  x  and  y  in 
Equation  A-l  results  in 

«  x    -  2        cxa+y2;/z      -     r     -      cos  e         3  A_2 


ry       z        (x^y^ya.    -  ~r   -   sin  e     ;  A_3 


a     -         Vx2-  _         jy_    _         sine 

Wx  ~       l  +(v/x)z  r2-  r  >  a-4 


a     _  .  Vx  _       x         _        cos  e 

WV  "       \  +  (Y/x)z  r2-  r  '   '  a- 5 


Now  consider  a  function    ^  Cr^ZL,©)        ,  wherein  r  and   0         are 
functions  of  x  and  y  through  Equation  A-l.    Taking  the  first  partial 
derivative  of      cf>    with  respect  to  x ,  y,  and  z  utilizing  Equations  A- 2 
through  A-5  gives 

<fiy.    -    ^r  rx  -V   0e  6*  =  ^r    C°S©  "  ^e    ^p^    ,  a-6 


</>y   -    $r    Py+     0o  By    -    ^r  5m0t    0 


cose 


r        >  a-7 


0Z    =      ^z    Zz     =     p'z         .  A- 8 
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The  second  partial  derivations  are  now  obtained  as  follows  through 

Equations  A-6  through  A-8: 

^             r    s                             ^         Sine    %   r  _/        ^^^^        ^        Sin  6    n 
^xx    =  C  ar    Cose  -    §^    — p — J  (  0r    COS6  -  fte  —  J 


A- 9 
S'm  e  cose    ,      0   ^  sine   cose 


-  0rr  Cose  +    0r    —    -v-  }2ee    p=r    -  Z  jdr  e 


2  0, 

fay^Ci?  Sine*  |^    ^^-)(  ^Sine  +  ^o   -^P-) 

cos2e    ,         /  cos2e 


0*-*-    5m20    +    0r      CO;    U    -v       0 


ee 


r  ^°"  r 

cose  sine  0   _/  cose    s'ms 


A- 10 


+  2^ro     L"bV  -    2^ 


9 


^zz      -       li^    (0O     -      0zz 


A- 11 


Adding  Equations  A- 9  ,  A-10  and  A-ll,  the  Laplacian  in  cylindrical 
coordinates  becomes 

V2    0     =      ^r    -V    -L-    ^r    +  "f*    0ee     +     ^z  A-12 


The  partial  derivatives  are  replaced  by  their  respective  central  difference 

2  2  ■        • 

operators     16     with  errors  of  order  h     where  h     is  the  square  of  the 

distance  between  adjacent  mesh  points  (See  Figure  A-l.)  resulting 


in 


A- 13 
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Figure  A-l 
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where  the  symbols  used  in  Equation  A- 13  are  shown  in  Figure  A-l. 

Salvadori  and  Baron  [16]  ,  as  well  as  any  other  good  text  on 
numerical  methods,  derives  finite  difference  operators  with  higher 
order  errors  ,  i.e.,  more  accurate  approximations . 

It  was  found  that  this  operator  cannot  be  used  for  a  mesh  point 
lying  on  the  longitudinal  axis  of  the  cylinder.    At  these  points 
Equation  A- 12  has  singularities  in  the  second  and  third  terms.     In 
this  study  circular  symmetry  was  assumed,  thus  Equation  A- 12  reduced 
to 

V     ^    --      ^rr     +y^r     t    ^7z         .  A- 14 

There  still  existed  a  singularity  in  the  second  term;  however,  the 
problem  under  investigation  also  possessed  symmetry  with  respect  to 
the  origin,  i.e .  ,  0r    =  0  at  r  =  0.     Smith  £9]    shows  that  by 

Maclaurin's  expansion 

2 


0rCr)     =     0r(o)+r    0rr    (o)+^ra0rrr   (0)+---      ,       A_x 


c 


however,  0r(O)  =  0.    Therefore,  as  r  approaches  zero  in  the  limit 

the  term      -~  (ft  r  in  Equation  A- 14  is  the  second  partial  derivative 

of       0        at  r  =  0  plus  higher  order  terms.     Equation  A- 14  can  be  written 
for  mesh  points  on  the  longitudinal  axis 


V2  0    =    2.  0rr    +      0zz 
or  in  finite  difference  operator  form 

+  -^2  \  0Cr,z+n  -  Z0(r/z)  +  0cr  >z:-Ol    • 


A- 16 


A- 17 
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APPENDIX  B 
DEFINITION  OF  COEFFICIENTS 

In  order  to  simplify  the  writing  of  Equations  11-35,  they  are  first 
expanded  in  order  to  utilize  the  sixteen  group  data  available  in  reference 
7    . 

A  typical  equation,  for  example,  group  three,  is  written  as  follows 


«_ii i 


A3(0 

— v 

A3<sO 


i        •>  =  ' ; 


-i^T-    •> 


3    2  ^^3 


A3(3) 


c 


3-^3     +    °~i^  >   3— *    A 


A3(4) 


_^ 
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V_^i ^ , 

c  A3<5^ 

CLJ2 v . 

A3Cig) 

This  equation  is  then  written  in  the  form 

[A3C4)  +  A3C3)  C4)]03(^^  •v-[A3(3)(-y^  +  jrh^l 


9^3  Cir  +  ljZ)  ■+  [  A3C3)  ("7?)]    03<r>Z-l)   +  •  •  •     =    AsCO  0,0,7) 

B-2 

+  A3  CZ)  02  Cr,z)  +  A3  (5)  03  (r,z)  4-   «  •   . 


+  A3 08)  0*0,2)  -v    S3Cr,z) 

The  general  form  of  this  equation  is 

[Ag  (9+1)  +  AS<g)  C^)]0gCr,z)  +  [  AgCg)  (f^  ~  )] 
0^O+\,z)  +  [  Ag  (<$)(   v>~  ETvTV]  0<3Cr-\,z)  +  [A3  (g) 

C  TT^l    ^3^r,Z+i)    4-    [  A<3Cg)(    y7*)]   ^gCr^Z-i)  B-3 

-  A3  CO    9$,Cr,z)  -t-  ...  -v  A3  (3-v)    ^g-^   Cr,-z^ 
-v  A3  (  3  +  2)  ^Cr.z)  +    -  -  •    +    A3  (18)  0,6  (r,Z) 
+   S3  (r,  z)    . 
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This  general  equation  is  condensed  further  to  the  form 

Al^gCr-ii-^)  +  A2  0gcrvz)  +  A3  09(r+i,z)+A40g(r,  z-0 
+  A5  09CT,  z+\)  =  AS,  0i(r}z")  +  •  •  •  +A-S^-i  0g-i  (r,-z.)     B_4 
+  AFg   ^9  Cr?z.:>  +  '  •  •    +  AF18   0,©  O,:z0  ■+  59Cr,zO  . 

In  energy  group  one,  Equation  B-3  does  not  apply.    The  reason  for  this 
is  that  the  terms  on  the  right  side  of  Equation  B-3  are  not  fully  applica- 
ble to  group  one.    These  terms  have  the  following  meaning:    the  sources 
of  all  neutrons  that  enter  group  g  due  to  down  scattering  from  all  higher 
energy  groups  plus  neutrons  born  into  group  g  due  to  fissions  in  all 
higher  energy  groups  are  represented  by  the  first  set  of  terms  on  the 
right  hand  side  of  Equation  B-3  |    /\o(  \)    ch  ,   (\~  "2.")    +-     •   •     • 

+  Ag  (9-1)  09-,   Cr,z)l  • 
In  group  one  these  terms  must  be  zero  since  there  are  no  higher  energy 

groups . 

These  apparent  inconsistencies  were  accounted  for  in  the  compute: 

program,  i.e.  ,  these  terms  were  set  to  zero  in  the  first  group  equation. 
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APPENDIX  C 
FINITE  DIFFERENCE  MESH  POINT  DESIGNATION  SYSTEM 

In  order  to  take  advantage  of  symmetry  in  the  study  and  hence 
reduce  the  computer  core  storage  space  needed,  the  following  mesh 
point  scheme  was  adopted:    a)  all  points       (fi  (l,x)  =         (ft  (3,x)  and 
d>  (x,  1)  =        d>  (x,3)  due  to  symmetry,  b)  all  extrapolated  boundary 
points         0  (L,x)  =         0  (M,x)  =  0  and  c)  no  mesh  points  have  zero 
index  due  to  computer  language  incompatability .    See  Figure  C- 1  for 
further  details . 
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Figure  C-l 
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APPENDIX  D 
GAUSS-SEIDEL  ITERATIVE  METHOD 

In  this  technique,  which  is  also  called  the  unextrapolated 
Liebmann  method,  the  point  single  step  iterative  method,   or  the  method 
of  successive  displacements,  the  latest  iterative  values  are  used  as 
soon  as  they  become  available. 

Consider  the  solution  of  Poisson's  equation 


V20     -   £(X,Y)     , 


D-l 


in  a  rectangular  region  as  in  Figure  D-l.    Assume  the  solution  is  proceed- 
ing column  by  column  from  left  to  right. 
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Figure  D-l 
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D-2 


If  the  point  at  which  the  equation  is  presently  being  solved  is 
<p  (i,j),  the  iteration  formula  is 


where  the  standard  five-point  central  difference  approximation  is 
substituted  for  the        V        operator. 


POINT  SUCCESSIVE  OVERRELAXATION 

This  method,  also  called  the  extrapolated  Leibmann  method,  can 
be  thought  of  as  the  sum  of  the  Gauss-Seidel  iterate  Equation  D-2  and 
the  nth  iterate.     It  is  written  as 

,  n+i  {-  l    r      "n  +  i  yi 

+    0^+V  »>  j-0  +  9^Ci,  i+O  +  h2SCi^]3  D"3 

where  w  is  the  overrelaxation  or  acceleration  factor. 

The  optimum  value  of  this  factor  in  general  is       I  <  wopt  <  2.       and 
can  be  computed  from 

W        =  ,  *- D-4 


opt  \-\  a/  i  -   x,s 


where        A,    is  the  largest  eigenvalue  of  the  problem  resulting  from  a 
solution  using  the  method  of  successive  displacements,  i.e.  ,  the 
Jacobian  method.    In  practice,  however,  the  eigenvalues  are  seldom 
known  and  finding  them  by  solving  the  problem  by  successive  displace- 
ments can  be  very  difficult  due  to  slow  convergence.    References     8,  9, 

10,    17     recommend  an  experimental  process  for  determining  w      ,  for  a 

opt 

particular  problem.     In  this  process  various  values  of  w  are  tried  until 
one  is  found  that  gives  the  most  rapid  convergence.    This  is  then  c 
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to  be  w       .    Reference     16     suggests  two  other  methods  for  determining 

w       . 
opt 

The  experimental  approach  described  above  was  used  in  this 
study. 

The  rapid  convergence  experienced  with  this  method  on  many 
different  problems  has  proved  far  superior  to  the  ordinary  Gauss-Seidel 
or  successive  displacements  method  without  overrelaxation.     References 
8,    17    contain  some  examples. 

References     8,  9,   10,    17     relate  more  explicit  details  on  the 
method  of  successive  overrelaxation  in  varying  degrees. 
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APPENDIX  E 
INPUT  DATA  SYMBOLS 


SYMBOL 

NAME  &  DESCRIPTION 

UNITS 

REMARKS 

GS 

MESH  SPACING 

CM. 

VALUES  FROM   .75 
TO  2  .  0  HAVE  BEEN 
USED  SUCCESS- 
FULLY 

GSSQ 

VALUE  OF  MESH 
SPACING  SQUARED 

CM.2 

C 

TOTAL  NUMBER  OF 
ELEMENTS  OR 
ISOTOPES  IN 
THE  ASSEMBLY 

UP  TO  FIVE  HAVE 
BEEN  USED 

R 

NUMBER  OF  MESH 
POINTS  IN  THE 
RADIAL  DIRECTION 
PLUS  ONE 

ADDITIONAL  POINT 
LIES  OUTSIDE 
REGION  OF  SOLU- 
TION BUT  IS 
NEEDED  FOR  DERIV- 
ATIVE BOUNDARY 
CONDITION 

Z 

NUMBER  OF  MESH 

POINTS  IN  THE 

AXIAL  DIRECTION 
PLUS  ONE 

IS  IN  EFFECT  ONLY 

HALF  THE 
AXIAL  DIMENSION 

SS(I) 

ENERGY 

DEPENDENT 

SOURCE 

NEU./CM." 

SPACE  DEPENDENCE 
MUST  BE  ACCOUNT- 
ED FOR  IN 
DIFFUSION 
EQUATIONS 

S(I) 

GROUP  RELATIVE 
MAGNITUDE  FACTOR 

DELE(I) 

ENERGY  GROUP 
WIDTH 

Mev 

H(I) 

ATOMIC  ABUNDENCE 
OF  EACH  ELEMENT 
OR  ISOTOPE 

ATOM/ 
CM.3 

SIGIN 
(U,K) 

INELASTIC  SCATTER- 
ING CROSS  SECTION 
OF  THE  Ith  ELEMENT 
FROM  GP.J  TO  GP.K 

barns 

■ 
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SYMBOL 

NAME  &  DESCRIPTION 

UNITS 

REMARKS 

SIGTR(IJ) 

THE  TRANSPORT  CROSS 
SECTION  OF  THE  Ith 
ELEMENT  OR  ISOTOPE 
AT  THE  Jth  ENERGY 
LEVEL 

barns 

SIGF(IJ) 

THE  FISSION  CROSS 
SECTION  OF  THE  Ith 
FISSIONABLE  ELEMENT 
OR  ISOTOPE  AT  THE 
Jth  ENERGY  LEVEL 

barns 

SIGC(IJ) 

THE  CAPTURE  CROSS 
SECTION  OF  THE  Ith 
ELEMENT  OR  ISOTOPE 
AT  THE  Jth  ENERGY 
LEVEL 

barns 

CHI  (I) 

NORMALIZED  FISSION 
SOURCE  FOR  THE  Ith 
ENERGY  LEVEL 

NU(IJ) 

AVERAGE  NUMBER  OF 
NEUTRONS  RELEASED 
PER  FISSION  OF  THE 
Ith  FISSIONABLE 
ELEMENT  OR  ISOTOPE 
AT  THE  Jth  ENERGY 
LEVEL 
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